OSU-ECE Report NASA 92-01 


Phase 2 

Semi-Annual Report on 


//j- a 3 


7'VJ./ 


NONLINEAR STABILITY AND CONTROL STUDY 
OF HIGHLY MANEUVERABLE 
HIGH PERFORMANCE AIRCRAFT 


(NASA GRANT NO. NAG-1-1081) 
Date: February 14, 1992 
R.R. Mohler, Principal Investigator 


Oregon State University 

Department of Electrical and Computer Engineering 
Corvallis, Oregon 97331-3211 

(503) 737-3470/3617 

(NASA-CR-189911 ) NONLINEAR STABILITY ANO N92-19841 

CONTROL STUDY OF HIGHLY MANEUVERABLE HIGH 
PERFORMANCE AIRCRAFT, PHASE 2 Semiannual 

Report (Oregon State Univ.) 64 p CSCL 01C Unclas 

G3/08 0071421 J 

Graduate Research Assistants: S. Cho, C. Koo, R. Zakrzewski 
Undergraduate Participants (NSF Support): D. Aaberg, J. Young 
Visiting Researchers: J. Dory*, J. Kurek**, A. Yagen* 


ADA- Israel Support 
International Exchange Board 



OSU-ECE Report NASA 92-01 


Phase 2 

Semi-Annual Report on 


NONLINEAR STABILITY AND CONTROL STUDY 
OF HIGHLY MANEUVERABLE 
HIGH PERFORMANCE AIRCRAFT 


(NASA GRANT NO. NAG-1-1081) 
Date: February 14, 1992 
R.R. Mohler, Principal Investigator 


Oregon State University 

Department of Electrical and Computer Engineering 
Corvallis, Oregon 97331-3211 

(503) 737-3470/3617 


Graduate Research Assistants: S. Cho, C. Koo, R. Zakrzewski 
Undergraduate Participants (NSF Support): D. Aaberg, J. Young 
Visiting Researchers: J. Dory , J. Kurek , A. Yagen 


*ADA-Israel Support 
* International Exchange Board 


TABLE OF CONTENTS 


Page 

1. OVERVIEW 1 

2. STATUS OF THE MODEL ALGORITHMIC CONTROLLER 6 

3. TIME-OPTIMAL CONTROL 8 

3.1 Introduction 8 

3.2 Switching-Time-Variation Method 10 

3.3 Approximation for Systems Affine in Control 13 

3.4 Computer Implementation 13 

3.5 Simulations 15 

3.6 Conclusions : 16 

4. REFERENCES 39 

APPENDICES 

A. List of Project Publications 

B. On Nonlinear Model Algorithm Controller Design 

C. Analysis of Nonlinear Stability Using Robust Stability Analysis for Linear Systems 


1. OVERVIEW 


This research should lead to the development of new nonlinear methodologies for the adaptive 
control and stability analysis of high angle-of-attack aircraft such as the FI 8 (HARV). The present 
progress report reviews project research over the first half of the second year. 

The emphasis has been on nonlinear adaptive control, but associated model development, system 
identification, stability analysis and simulation is performed in some detail as well. Table 1 summarizes 
various models under investigation for different purposes. 

Models and simulations for the longitudinal dynamics have been developed for all types except 6 

in Table 1. A very preliminary analysis has been made on type 6 (neural net models) for adaptive control 

( 

thus far. It has been shown that dynamic accuracy roughly increases with ascending order of model type 
from 1 to 7, except that perhaps 3 (Volterra series) and 6 (neural nets) should be interchanged. However, 
such comparisons depend on how the models are utilized. Here, the focus is on adaptive control, 
generated by model-reference types 1 to 6, of a complex nonlinear aircraft motion represented by 7 
(nonlinear ordinary differential equations). Preliminary analyses use a nonlinear second-order 
approximation [1] which we found useful for changes in angle of attack (ex) by about 10°. A fifth-order 
nonlinear longitudinal model with the traditional stability derivatives generated as functions primarily of 
a, for a given altitude and mach number, successfully mimicked F18 flight trajectories [2], and is being 
utilized for our nonlinear adaptive-control studies at the present time. These models are discussed in the 
project’s first annual report [3]. 

Briefly, studies completed indicate that nonlinear adaptive control can outperform linear adaptive 
control for rapid maneuvers with large changes in a. Figures 1 and 2 compare the transient responses 
where the desired a varies from 5° to 60° to 30° and back to 5° all in about 16 sec. Here, the 
horizontal stabilator is the only control used with an assumed first-order linear actuator with a 1/30 sec 
time constant. Unfortunately, an additional rate constraint significantly reduces the system performance 
for both the nonlinear and linear adaptive control as shown in Figures 3 to 5 and analyzed in the next 
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Table 1 . Aircraft Models 


Type 

Purpose 

Remarks/Limitations 

1. Linear perturbations at 

a = 5°, 15°, 35°, 60° 

Local control, check of nonlinear 
system, application of well devel- 
oped linear control methodologies 

Local stability 

Only valid for small maneuvers 
Special case of types 2-5 

2. Gain scheduled (non- 

linear function of a) 
from 1 

Gain-scheduled adaptive control 
based on well developed methodolo- 
gies 

Simplified description of complex 
system 

Approximate stability 

May have stability problems 
with small number of reference 
states and/or large fast maneu- 
vers 

3. Volterra series 

a) at reference states 

b) general case 

Nonlinear adaptive control via cross- 
correlation and/or i priori dynamic 
structure 

stability approximation 

Simplified dynamic description of 
complex system 

Non-orthogonal series approxi- 
mation 

Sufficiency of 2 or 3 kernels 

Large computation time for 
adaptation 

4. Bilinear system 

a) continuous 

b) BARMA 

5. Polynomial time series 

Nonlinear adaptive control via model 
reference identification (NLMRAC) 

Stability approximation 

Simplified dynamic description 

Large computation time 

Bilinearizing controllers may be 
more practical than linearizing 
ones 

Polynomial approximation may 
be more accurate but more time 
consuming than linear or bi- 
linear approximation 

6. Neural network 

Potential application to adaptive 
control 

Probably less accurate than 4 or 
5 for a given data set but accu- 
racy may be more robust out- 
side the available data set 

7. Nonlinear ordinary 

differential model 

Accurate approximation to fast large 
maneuvers for 'final* design and 
simulation 

Stability 

Neglects flexible modes and 
other complications 
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Figure 4. Response of Nonlinear MAC with Stabilator Limit of 40° per Second 
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Figure 5. Response for Nonlinear MAC with Stabilator Limit of 40° per Second and Slower Reference 
Trajectory 

section. Such lack of controllability can be improved, of course, by introducing thrust vectoring as used 
in [2]. Appropriate thrust-vector control to supplement the traditional pitch-motion stabilators is 
underway for the nonlinear adaptive controllers, and preliminary results are encouraging. 

A preliminary analysis of time-optimal control of a is studied in Section 3. Here, a new algorithm 
is derived from the switching-time variational method [4,5] and then applied successfully to the simplified 
second-order nonlinear model [1]. The method is presently being adapted to the more complex nonlinear 
fifth-order model. This study should provide a "yard stick" by which to evaluate controller performance 
as well as provide a base for more effective controller designs. As a byproduct of this analysis the 
complicated Jacobian of the longitudinal dynamics will be computed as a function of a and other 
variables. While it is used here to compute bang-bang controller switching times, it may have other uses 
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for approximate dynamic-system identification beyond the usual time-invariant linearized models at trim 


states. 

Finally, linear stability arguments are developed in Appendix C which tend to at least an 
approximation of the admissible range of model parameters as applied to the nonlinear second-order 
approximation [1]. 

2. NONLINEAR MAC ALGORITHM STATUS 

Model algorithmic control (MAC), described in [3], starts with 

+ + («0O - “.nodOO) M 

where 

« P. T 4>(k) ( 2 ) 

4>(k) * [a,a 2 ,a 3 ,q,qa,qa 2 ,qa 3 ,u,ua,ua 2 ,ucc 3 ,l] T (k) 

As the control at the moment k must be already computed at moment k the values of a(k) and q(k) are 
not available for its computation so their estimates must be used instead. The correction term is taken 
to be the prediction error from the moment k-1 and the equation becomes 

a.rfOc+l) = “.nodfr+l) + («(k-D - “^(k-D) W 

with 

®mod(k + l) “Pjfa) 

|>(k) = [d,d 2 ,d 3 ,q,qd,qd 2 ,qd 3 ,u,ud,ud 2 ,ud 3 ,l] T (k) 
a(k) = p a T 4>(k-l) + (a(k-l) - a^Ck-l)) 
q(k) = p^Ck-l) + (q(k-l) - q^Ck-l)) 
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The controller is assumed to know the values of angle of attack and of pitch rate at the moment k-1. 
Then it estimates their current values a(k) and q(k) taking into consideration previous prediction errors, 
and based on them it calculates the control required to achieve o^f at the moment k+ 1. The value of 
control is found as: 


u(k) = - Pla< * ~ p2 ‘ &2 ~ p3g&3 ~ p4a ^ " p3a ^ & ~ P<5g ^ &2 ~ P? « 4&3 ~ Pl2a 

Pg« * P 9 .«Pio -“ 2 + Pu .« 3 


( 5 ) 


where 


& t = *«&+» - (“Oc-D - cc^Ck-l)) (6) 

and a = a(k), Q = Q(k) as described above. 

This algorithm was made to be adaptive, or self-tuning, by incorporating on-line identification of 
the parameters. A recursive least squares (RLS) algorithm was implemented in the form taken from 


P(k) 


Q(k-2) <t>(k-l) 

X(k-1) h. <j>(k-l) T Q(k-2) <J>(k-l) 


( 7 ) 


Q(k-1) 


1 

X(k-1) 


Q(k-2) 

\ 


Q(k-2) <t>(k-l) <{>(k-l) T Q(k-2) ' 
A(k-1) + <j)(k-l) T Q(k-2) <j>(k-l) ; 


( 8 ) 


e(k-l) = y(k) - p T 4>(k-l) ( 9 ) 

where y may denote a or q and p may stand for p a or p q , respectively. The forgetting factor X was 
introduced to enable the algorithm to change the estimates of parameters with the change of operating 
conditions. To avoid the unlimited growth of covariance matrix Q at the steady state when the input is 
not persistently exciting, the variable forgetting factor policy was implemented: 

X(k) = 1 - e ^ (10) 

e(k) 2 
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where e(k) is the current prediction error, e(k) is the average prediction error from last 10 samples, and 
e is equal to 0.01. As an additional precaution, the trace of the covariance matrix Q was monitored and 
Q was reset to diagonal matrix whenever the threshold value was exceeded. 

To further damp the response, the controller is designed to minimize the one step ahead cost 
function: 

J » (y^Gc+D - y^i}) 2 + p m - u(k-i)) 2 (ID 

with y mod , y r as before. Minimization of (11) with respect to u(k) yields 

m , : a ) b * PU(1I - 1) ' (i2) 

b 2 + p 

where 

a • Pta« + P 2a “ 2 + P3«<* 3 + M + Ps«q« + P6«<I« 2 + P7.q« 3 + Pl2a 
b « P 8a + P 9 « a + Pl0«* 2 + Pll«* 3 

Obviously, for p = 0 (12) reduces to (5) while for p = oo we have u(k) = u(k-l) = const. 

This controller is used in Figures 2, 4, and 5 with only the linear portion of a,b used in Figures 
1 and 3. The algorithm will be generalized to include thrust vector control and variable-horizon cost. 

3. TIME-OPTIMAL CONTROL 
3.1 Introduction 

Various control strategies have been developed by the team and to find their merits it seems usefu 
to have an idea of what are the best output and state trajectories theoretically possible, given the existing 
constraints on the control variables. For substantially nonlinear systems the problem of synthesis of the 
optimal feedback control law is usually untractable. On the other hand, there exist numerical technique: 
that allow us to calculate "open loop controls" - i.e., the specific control signals necessary to achieve the 
minimum performance index. Aware of the difficulties connected with the controller synthesis problen 
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we do not seek its exact solution; at this time, we merely want to find the limit for the performance of 
a controller assuming perfect knowledge of plant dynamics and absence of any unforeseen disturbances. 

This report is concerned with the problem of time-optimal control in which we are interested in 
transferring the system’s state from an initial value to some prescribed terminal set in minimal time. In 
the aircraft problem this might mean changing the flight’s pitch angle, path angle, or angle of attack from 
an initial equilibrium value to some other terminal value, preferably also with all other states moving to 
the equilibrium. The control value (stabilator or elevator angle) is naturally bounded from below and 
from above. For some systems it turns out that in case of such simple cube-type constraints on control 
variables, the time-optimal control is of bang-bang type. However, for quite a large class of systems that 
are affine in control, we may approximate any measurable control signal with a bang-bang signal with 
arbitrary accuracy in the sense that corresponding state trajectories are arbitrary close to each other in 
L 1 metric. Hence, also time-optimal control, if it exists, may be approximated by a bang-bang control, 
even if it contains singular arcs. Therefore, the approach presented here is to find the bang-bang control 
that will minimize the transition time. The computational algorithm used here is the switching-time- 
variation method developed in [4,5]. Since the algorithm gives as an output a control signal with finite 
number of switchings, it is tacitly assumed that with large enough finite number of switchings, we are 
able to achieve good enough approximation of optimal control. This, unfortunately, does not follow from 
the theory I am aware of, since the above mentioned approximation result holds only for bang-bang 
signals with possibly infinite or even uncountable number of switchings. This delicate question is left 
aside for the time being to be clarified later. Another point worth indicating here is that resulting control, 
in an attempt to approximate a continuous "singular" control, may have inter-switching times very small, 
thus precluding any practicality of the approximation. This, however, is of no concern to us since, as 
mentioned before, we are interested only in finding the best possible output, or state, trajectories - not 
the actual control signals corresponding to them at this time. 


9 


A computer program has been developed for numerical solution of the problem. The program, due 
to its modular construction, easily allows various plant models to be plugged into it. The switching-time- 
variation method is used in it for fixed terminal time with the quality function being the weighted distance 
of the target set. Then the smallest such time is found that allows it to hit the target exactly, and finally 
the optimal number of switchings is iteratively found that gives minimal transition time. 

In what follows the switching-time-variation method is briefly characterized in Section 3.2. Section 
3.3 discusses briefly the approximation theorem for bang-bang controls in systems affine in control. The 
computer implementation of the algorithm is discussed in Section 3.4. Section 3.5 contains the test 
results of the program for a second-order model of longitudinal dynamics of an aircraft. The concluding 
remarks discuss the possibilities of application of the computer package to solutions of more complex and 
problems more close to reality. 

3.2 Switching-Time-Variation Method 

The switching-time-variation method used here was taken from [4], and the original thesis [5] was 
also consulted for the details. The method is designed for the computation of optimal control in the class 
of bang-bang control signals with finite number of switchings. The quality criterion is assumed to be 

J = Q (f„0O + go(x)u(t))dt (13) 

for the system of the form 

^ = f(x) + g(x)u(t) (14) 

dt 

where x e R n , u e R 1 , t e (tg.tf]. To ensure the existence and uniqueness of solutions of (14), f and a 
are assumed to be continuously differentiable with respect to x. The control values are constrained by I 

-1 s u(t) si (15) I 
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Of course, any control constraints of the cube-like type u^,, ^ u(t) ^ u^^ may be transformed to form 
(15) for system affine in control. The control objective is to minimize the quality criterion (13) for given 
initial state xq with possible penalty term connected with final state already included in f 0 and go by 
standard transformations, assuming that admissible controls are bang-bang with finite number of 
switchings. The version of the algorithm described in [4] was developed for systems with scalar controls 
and the computer program described here is also designed for this special case. However, it is not of 
any particular difficulty to generalize the algorithm to the case of u e R m . If the need arises, the 
computer program may also be modified to accommodate this possibility. Here the scalar version will 
be presented because of its notational simplicity. 

The method is an iterative one - in each step the gradient of the quality criterion with respect to 
switching times is computed. The switching vector is defined as: 

x = (t (16) 

where N is the number of switchings, with constraints: 

1^ s t, s - s t H s t, (17) 

The control value on the interval [r ; ,r i+1 ) is then equal (-1) 1 . The augmented system is defined as: 


^ = f(x) + g(x)u(t) 
at 

where F = (f 0 ,f T ), F = (g 0 ,g T ), Xq = (0,Xg), and the adjoint system equation is 


dX 

dt 


,|<0 + |f(t)u(t) 

^ dx uX. j 


(18) 


(19) 


with terminal conditions X^tf) = (dJ/3jq)(tf). Then the gradient of the quality criterion with respect to 
the switching vector may be calculated by means of the formula 
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( 20 ) 


H- - 

with function <£ defined by 

4>(t) = 2(g(x(t)),A.) (2D 

with gradient calculated the method consists of iterative descent steps 

r.(k+l) = ^(k) (22) 

dti 

where lq are such that constraints (17) are satisfied and sufficiently small to ensure that J(k+ 1) < J(k). ’ 
The algorithm is terminated if either the gradient is zero or no feasible (i.e., descent) step may be 
executed. 

On top of the algorithm of finding the optimal switchings with their number given there is an outer 
loop modifying this number. If the optimal control results in a constraints t ; <. r ;+1 active than the 
switchings i and i+ 1 should be removed. On the other hand, if there are two zeros of 4>( t) not coinciding 
with any of the switching times than two switchings should be added between these zeros. After the 
modifications of the dimensionality of the switching vector the inner loop of optimization is again 
performed and the process is terminated when no more changes of the number of switches are necessary. 

It is worth noticing that the above algorithm of finding the optimal bang-bang control may be also 
generalized for broader class of systems dx/dt = f(x(t),u(t)). The main difference would be the formula 
for function <f>(t0) which would become 

<j>(t) = ((f^tO.^W) - ^(O.u^)) , X) (23) 

Of course, the technical assumptions ensuring the existence of solutions should be satisfied. 
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3.3 Approximation for Systems Affine in Control 

The algorithm described above calculates the optimal control within the class of bang-bang control. 
However, for systems affine in control a result is available stating that we may approximate an arbitrary 
admissible control with a bang-bang control such that corresponding trajectories are arbitrary close. 

The theorem, stated and proven in [6], assumes that we have a system of the form (14) with 
constraints (15). Functions f and g are continuously differentiable, and a Lipshitz type condition 
< f(x) + g(x)u,x> < K(1 + (| x | 2 ) preventing finite escape-time is also assumed to be satisfied for all 
x in the region of interest. Then an arbitrary measurable control signal u(t), t e [tg,tf] satisfying (15) is 
considered with corresponding state trajectory x(t). Then the theorem states that given any e > 0 is 
always possible to find a bang-bang control u*(t) satisfying | u*(t) | ^ 1, such that the corresponding state 
trajectory x*(t) approximates x(t) uniformly on [tg,tf] with accuracy less than e, i.e., |x(t) - x*(t)| <, e 
for all t e [tQ.tf]. 

Although the theorem stated above considers a bang-bang control with not necessarily finite or even 
countable number of switchings, it gives some justification to using the switching-time-variation method 
for systems with singular optimal controls. Intuitively for reasonably smooth systems there should be 
some kind of continuity enabling in turn approximating the bang-bang control u* with a sequence of bang- 
bang signals u with finite number of switchings. However, I am not aware of any such result, and in 
monograph [7] from 1990 the aforementioned result is cited after [6] as the only available. It still seems 
feasible to come up with some, maybe more restrictive, assumptions which would justify using finite 
number of switchings. 

3.4 Computer Implementation 

The algorithm discussed in Section 3.2 was implemented as a quite general software package. It 
finds the time-optimal control for the case when the terminal set is a single point y. The time-optimal 
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problem with fixed terminal set is replaced with a sequence of fixed time and free terminal state problems 
with quality index 


i - 2 p««,) - yf P 4 > 

Switching-time-variation method is used to solve this problem, and the desired final time is decreased if 
the resulting quality is zero or is increased in the opposite case. This iteration is repeated until we get 
to the limit time t* below which the quality is always positive, i.e., it is not possible to find a bang-bang 
control transferring the system from Xq to y. 

The optimization method described in Section 3.2 was modified somewhat in details of the gradient 
minimization routine. Instead of performing single step in the direction, a directional search is performed 
with constrained step size. A combination of two-point gradient parabolic approximation and three-point 
non-gradient parabolic approximation is used to find the minimum in the direction. The generation of 
the descent direction is also somewhat different. First, if any of the constraints (15) are active and the 
gradient points outside the feasible region, the gradient is projected on the proper constraining 
hyperplane. The special structure of constraints causes the projection to consist solely of putting the 
appropriate coordinates of the gradient to zero. Then the direction is tangent to the constraining 
hyperplane, and we get an optimization problem of reduced dimensionality. This problem is solved using 
a conjugate gradient method in the version proposed in [6], The conjugate gradient is restarted not only 
every N iterations, where N is the current dimensionality of the problem, but also whenever the set oi 
active constraints changes - i.e., when the algorithm hits or leaves a constraining hyperplane. The 
termination of the procedure occurs when the projected gradient is zero - i.e., no feasible descent step 
is possible, or equivalently when the dimension of the current optimization hyperplane becomes zero. 

The calculation of the quality criterion and of its gradient involves numerical integration of Eqs 
(18) and (19). This is done using a fourth-order Runge-Kutta integration method. To integrate th< 
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adjoint equation (19) the whole state trajectory resulting from integrating (18) must be stored, but for 
calculation only a small number of points from the costate trajectory is needed. 

The program is written in the fashion enabling easy substitutions of different plant models and 
different optimization tasks. To use another model one has simply to provide the routines calculating the 
right-hand sides of Eqs. (18) and (19). The problem is defined in a straightforward fashion by setting 
the values of initial state, terminal state, initial estimate of the final time, etc., in the main routine. The 
whole program is written in C programming language, and although compiled and run on an IBM PC, 
it may be easily ported to any machine with C language compiler. The only difficulty that may occur 
with more complex systems is the rather severe storage requirements - whole state trajectory has to be 
stored with sufficiently small discretization step in order to calculate the gradient. And, of course, there 
will always be the problem with the speed of calculations for higher dimensional systems. 

3.5 Simulations 


The program described above was tested on a model previously used (in our NASA project), i.e., 
a simplified longitudinal aircraft model of second order (so called "Stalford model") described in [6]. 

For the aircraft model the problem solved was to increase or decrease the value of angle of attack 
with the requirement that the maneuver should move the system from the equilibrium corresponding to 
starting value of the angle of attack to the equilibrium corresponding to its final value. The control 
signal, the elevator angle, was assumed to be between 0 and -20 degrees. The series of maneuvers 
simulated was labeled in the following way: 


maneuver A: 
maneuver B: 
maneuver C: 
maneuver D: 
maneuver E: 


from 0 to 15 degrees; A’: from 15 to 0 degrees 

from 0 to 18 degrees; B’: from 18 to 0 degrees 

/ 

from 0 to 20 degrees; C’: from 20 to 0 degrees 
from 5 to 15 degrees, D’: from 15 to 5 degrees 
from 5 to 18 degrees, E’: from 18 to 5 degrees 
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maneuver F: from 5 to 20 degrees, F’: from 20 to 5 degrees 
maneuver G: from 10 to 15 degrees, G’: from 15 to 10 degrees 
maneuver H: from 10 to 18 degrees, H’: from 18 to 10 degrees 
maneuver I: from 10 to 20 degrees, I’: from 20 to 10 degrees 
maneuver J: from 15 to 18 degrees, J’: from 18 to 15 degrees 
maneuver K: from 15 to 20 degrees, K’: from 20 to 15 degrees 

The results of the optimization for each of these maneuvers is depicted in Figures 6-16. It may be 
observed that for all of them the time-optimal control had only one switch. In all cases the time-optimal 
trajectory for a had a substantial overshoot and consisted of an almost linear first portion with high slope 
before the switch and of slowly decreasing second portion. 

3.6 Conclusions 

The computer program presented here is suitable for calculating the time-optimal controls for 
arbitrary finite-dimensional systems which are affine in control. The simulation results discussed here 
have mainly testing value showing that the program is in operation. The next step should be to use the 
program on some more complex models such as the fourth-order longitudinal-aircraft model, which 
together with a linear actuator is definitely affine in control. The resulting time-optimal trajectories for 
different maneuvers could be used as benchmark tests for other controllers or as reference trajectories 
for time-series-based, adaptive, one-step-ahead (or many-steps-ahead) control. The work on this is 


underway. 
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ON NONLINEAR MODEL ALGORITHMIC CONTROLLER DESIGN 
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1. INTRODUCTION 

Two nonlinear algorithmic controllers, MAC, are studied here. One uses a 
block-canceling Volterra approximation, and the other MAC consists of solving an 
approximating polynomial time series instate and control. Both methods synthesize 
discrete control sequences and are applied successfully to the control of a simple 
nonlinear longitudinal aircraft model for large variations in angle of attack. 

The Volterra-series approach used here was introduced by Modyaev and Averina 
[1], and a form of inverse generating control according to an assumed structure is 
presented by Harris [2]. This work formed the basis for the methods used here. The 
high angle-of-attack aircraft model derived by Stalford, et al. [3] was the plant 
simulated for the MAC application. In many traditional design studies, a sequence 
of linearized perturbation models are derived for different equilibrium flight 
conditions with linear controllers appropriately derived. Linear adaptive control can 
be derived according to nonlinear gain scheduling of the control law. A highly 
successful version of such control, which includes proportional plus integral plus filter 
(PIF) terms, is presented by Ostroff [4,5]. However, such designs usually require a 
large number of set-point design computations, and may have stability problems for 
large fast changes in angle of attack and/or mach number. 

For generation of the nonlinear control, a nonlinear time-series based model 
reference is used. In order to identify such model, experimental data was collected 
for angle of attack (a) and pitch rate (q) subject to random steps of control 


(stabilator, 5). To capture such phenomena as limit cycles in the data the steps were 
rather long (40 s). There were 64 such steps with time discretization of 0.1 s 
resulting in 25,600 points in a state plane for 64 values of control. 

For a least-squares simulated data Fit, the following approximation was 
surprisingly accurate: 

®(k+1) = Pi*®(k) + P 2a « 2 (k) + P 3a « 3 (k) + 

P^k) + P 5a q(k)a(k) + Pe a q(k)a 2 (k) + P 7a P(k)a 3 (k) + 

P 8a u(k) + Pg a u(k)o(k) + p 10a u(k)a 2 (k) + p 11a u(k)a 3 (k) + p 12a 
q(k+1) = P 1q a(k) + P 2 q® 2 (k) + P3q® 3 (k) + 

P4qQ(k) + P 5 qq(k)a(k) + Pe<,q(k)a 2 (k) + p 7q q(k)a 3 (k) + 

P*,qu(k) ♦ PgqUfkJaW + p 10q u(k)a 2 (k) + p 11q u(k)a 3 (k) + p 12q 

Even limit cycles are accurately rendered by this model, as well as the stable zone 
behavior, although large discrepancies occur when the control values are close to the 
stable/unstable zones border. 

2. ADAPTIVE CONTROL APPROACHES 
2.1 Nonlinear Volterra-Based Control 

Here, as in [6], the Volterra series serves as a conceptual starting point for a 
nonlinear time series base control. Continuous time controllers based on Volterra 
series were systematically developed in [7] with formulae for the controller’s kernels 
given those of the plant and of the desired feedback system. In particular, the 
problem of so-called exact feedback linearization was solved here. However, those 
formulae may be of limited practical value because of the properties of Volterra 
series under feedback. The problem is that even finite (e.g., second order) Volterra 
series of the open loop results in infinite Volterra series of the closed loop. This 
makes it necessary for the controller to include theoretically an infinite number of 


compensating terms even for a quadratic system. The same problem for the discrete 
time systems was treated in [1] with multidimensional Z transforms to derive the set 
of formulae equivalent to those for so-called exact feedback linearization [8]. 
However, they also provided a very elegant transformation of which results in a 
controller requiring only as many Volterra terms as there are in the assumed plant. 

One attractive feature of this controller is that its structure makes it possible to 
utilize it not only with models represented in the form of Volterra series, but in fact 
with any model with easily divided linear and nonlinear parts of the dynamic 
equations such as (2) above. 

The following algorithm results: 

a) according to the linear part of the plant, calculate the linear control u L (k) 

b) calculate the predicted value of the output at the moment k 

y(k) = L (y(k- 1 ),...,y(k- M),u(k- 1 j u(k-M)) 

Af(y(k-1) y(k-M),u(k-1) u(k-M) 

c) solve the "linearizing" control equation for x(k) such that 

W(y(k).y(k-1) y(k- M + 1 ),u L (k) -x(k),u(k- 1 ) u(k-M+1)) = 

=Z.(x(k),x(k-1 ),...,x(k-M + 1 ),?(k),y(k-1 ),..., y(k- M + 1 )) 

3) calculate the control by 

u(k) = U L (k) - x(k) 

This algorithm becomes a sort of prediction controller which tries to estimate the 
effects of the previous controls knowing the previous values of outputs and then to 
adjust the current value of control so that the nonlinear part of predicted output is 
canceled. 

This discrete time nonlinear a control algorithm is generated according to an off- 
line identification of model (1) with a nonlinear aircraft simulation based on [3]. 
Also, a linear controller was designed according to the linear parts of (l)-(3). 

The design was performed to obtain the closed loop model reference behavior 
of the form 


-l 


G(z) = 0.05 /(z 2 - 1.6z + 0.65) 


In order not to cancel the zero of the plant, the observer polynomial (z-0.7) was 
introduced. The algorithm for the control value u(k) is as follows. First the estimate 
of the output at moment k is calculated from (1) with k replaced by k-1. 

Then it can be shown that the control becomes 

u(k ) + P8 « u lQ<) - (p 2a « 2 * P3«« 3 + Ps««q + Peqq« 2 + p 7 «q« 3 + P 12 .) (4) 

(P*« + Pg«“ + Pio«“ 2 + Pn«® 3 ) 

with a(k) and 3(k) designating estimates taken from (1). It is seen that if there are 
no nonlinearities in the model the control reduces to a regular linear controller 
u = u L . 

Simulations were run to test the controller performance especially in the 
unstable range of angle of attack. The system is successfully stabilized and the 
transients are very smooth and without significant overshoots for the nonlinear 
control as demonstrated by Figure la. By different choice of the reference model 
it is possible to obtain much faster, but at the same time much more "nervous” 
transients. The elevator control is also relatively smooth and within the range 
corresponding to the terminal equilibria. As can be seen from Figure lb, the similar 
linear control is unstable. 


2.2 On-Line Adaptive MAC Algorithm 

.1 

Model algorithmic control (MAC), described for example in [2], consists of 
solving the model equation for the value of control necessary to obtain required 
value of output. Usually this desired output trajectory is generated form the setpoint 
by means of a reference model. In case this model is linear, the algorithm in essence 
becomes a linearizing one. 


Here, the controlled output is assumed to be the angle of attack such that the 
reference equation becomes: 

a ref( k+1 ) = ®mod( k+1 ) + (“(k- 1 ) - “modCk" 1 )) (5) 

with 

«mod(k + 1) = P« T ^(k) 

<M k ) = [a, a 2 , a 3 , q, qa, qa 2 , qa 3 , U, ua, Ua 2 , Ua 3 , l] T (k) 

) «(k) = P.V-1) + (a(k-1) - a^k-1)) 

q(k) = P q T 4>(k-1) + (q(k-1) - q mod (k-1)) 

The controller is assumed to know the values of angle of attack and of pitch rate at 
the moment k-1. Then it estimates their current values a(k) and q(k) taking into 
consideration previous prediction errors and based on them it calculates the control 

■-% 

required to achieve a ref at the moment k+1. The value of control is found as: 

U (k) = *' ~ Plg ° ~ P2 ““ 2 ~ P3a ° 3 P4g ^ ~ Psa ^* ~ P6a ^ £2 ~ P7 “^ 3 ' Pl2 ‘ 

Ps« + Ps«« + Pio«« 2 + Pii«“ 3 '* I 

< 6 > ■ 4 

where 

“r = «re((k + 1) - (a (k-1) - a mtxl (k-1)) 

and a = a(k), q = q(k) as described above. 

The results of the simulations are seen in Figures 2a, b. The-reference trajectory 
was chosen to be l/zM.6z+0.65). The actual output of the plant is seen to follow 
the reference very closely, even though the region of operation was that of the most 
severe nonlinearities. The control action is also remarkably smooth. 

The discrete time nonlinear state space model (1) describes the behavior of the 
complex nonlinear plant quite accurately in the entire region of operation. In 
practice, however, such a global model is rather difficult to fit, and consequently one 


should look for local approximations, depending on the current operating conditions. 
In such a situation, on-line adaptive control seems to offer an ideal solution. 

The algorithm discussed in the previous section can be made adaptive, or self- 
tuning, by incorporating on-line identification of the parameters. A recursive least 
squares (RLS) algorithm was implemented in the following form taken from [8]: 


P(k) 


Q(k-2) 4>(k-1) 

l(k-1) + <j>(k-1) T Q(k-2)4>(k-2)4>(k-1) 


e(k) 


(7) 


Qflc-1) = 1 f Q(k -2) - Q (k , r 2)^k- 1 ), <|>( k- 1^Q(k-2 X | w 

Mk-1) i A.(k-1) + 4>(k-1) T Q(k-2) <J>(k-1)J 


e(k-1) = y(k) - p T <J>(k-l) 


0 ) 


where y may denote a or q and p may stand for p„ or p q , respectively. The forgetting 
factor X was introduced to enable the algorithm to change the estimates of 
parameters with the change of operating conditions. To avoid the unlimited growth 
of covariance matrix Q at the steady state when the input is not persistently exciting 
the variable forgetting factor policy was implemented: 


A(k) = 1 


e (k) 2 


(10) 


where e(k) is the current prediction error, e(k) is the average prediction error form 
last 10 samples and e is equal to 0.01. As an additional precaution the trace of the 
covariance matrix Q was monitored and Q was reset to diagonal matrix whenever the 
threshold value was exceeded. Starting values of parameters were taken to be as in 
( 1 ). 

Figure 2 displays the simulation results for a reference model specified as 
l/(z 2 -1.8z+0.82). Remarkably exact following of the reference trajectory may be 
observed, although, surprisingly enough, the performance is slightly worse than in the 
nonadaptive case. Most probably this is due to the fact that prediction error now 
changes much more quickly because of the ongoing identification process. Thus, 
approximating the term (y(k+ l)-y tnod (k+ 1)) by (y(k-l)-y nx>1 (k-l)) may worsen the 


behavior of the system as two values of y ^ no longer correspond to the same 
parameter vector. Since, the on-line identification process assures (at least in 
principle) that the prediction error should asymptotically converge to zero it is 
possible that the correction terms in a(k), q(k), and in control equation (5) ought to 
be omitted. 

The performance of the adaptive nonlinear MAC controller was compared to 
the linear one, which uses the same control strategy but with a strictly linear model 
being identified and used for the calculation of the control action. Clear difference 
between the performance of linear and nonlinear controller can be seen from Figure 
3, particularly in control action at the setpoint a = 15°. The linear identifier has 
obvious difficulties with fitting the parameters of a linear model to the behavior of 
the plant which is highly nonlinear in this region. As a result, the control starts 
oscillating for a while. Also, it was seen that the nonlinear algorithm results in 
control plots that are more smooth, although they still contain one-pulse spikes. To 
eliminate these spikes weighting of the increments of control can be introduced into 
the algorithm with little performance deformation. 


4. CONCLUSIONS 

The nonlinear control applications to high angle-of-attack aircraft, as reported 
here, is of a preliminary nature. However, the analysis does suggest that nonlinear 
adaptive control can be quite effective to stabilize large rapid maneuvers in angle of 
attack. Of the comparisons made, the on-linear, nonlinear-time-series and adaptation 
performed the best and was quite superior to a similar linear MAC. 
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Figure 3a: Linear adaptive MAC (with Figure 3b: Linear adaptive MAC 
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Summary 

The stability analysis of an airplane using its nonlinear model is 
presented. The analysis is based on the robust stability analysis 
approach for linear systems. Then, based on analysis, a small 
static feedback gain is designed such that the robustness of the 
closed-loop nonlinear system stability is significantly improved. 
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1 . INTRODUCTION 


The stability is one of the most important issues in the cont 
system design. Recently there has been observed a great inter 
in the methodology of robust stability analysis and design 
robust control systems for linear dynamic systems [6] . 
objective of this paper is to investigate the applicability 
this approach for nonlinear dynamic system such as an aircr 
flight in high angle of attack/sideslip flight. The unsta 
control system can result in the plane crash. 


There is considered stability of nonlinear, simplified howev 
model of the airplane. The organization of the paper is 
follows. In section 2, the model of the plane is present 
Stability of the aircraft is considered in section 3. Final 
concluding remarks are given. 

2. THE AIRCRAFT MODEL. 

Model of an airplane is highly nonlinear, [4,5]. There 
usually, however, used simplified models for control sys 
design, e.g. [1,3,9]. In this paper we consider very simple mo 
given in [8] : 

x=A (x) x+Bu+D 


where x= 


is a state vector, a is the angle of attack 
degrees, ~q” is the pitch ratio in degrees per second and u is 
elevator control in degrees, 


9.168c (a) 1. 

-5. 73 o/ 


B= 


-1.83361 f-5. 473296' 

-8.5950J' u '[ 2 . 865000_ 


and c (a) is a nonlinear function. This function can, however, 
z 

approximated as follows: 


2 


-0.072815870 

for 

0° 

< 

a 

< 

14.74° 

0.088470922-2.3774/a 

for 

14.74° 

< 

a 

< 

17.40° 

0.033099050-1.4068/a 

for 

17.40° 

< 

a 

■< 

18.87° 

-0.016633734-0.4743/a 

for 

18.87° 

< 

a 

■< 

28.00° 


It is easy to find that 


-0.072818087 < c (a) 

z 

s -0.048161261 

for 

14.74° < 

a £ 17.40° 

-0.047751524 < c (a) 

z 

=s -0.041453149 

for 

17.40° < 

a s 18.87° 

-0.041768869 < c (a) 

z 

=s -0.033573019 

for 

18.87° < 

a s 28.00° 

This model approximates 

model taken 

from 

measured 

wind tunnel 


values of the T-2C airplane [7] . It is known that numerical values 

of c (a) and b 0 are uncertain, 
z 2 

Our purpose is to consider stability of system (1) in the range of 
angle of attack 0°ia^28°, and to find a static feedback which can, 
eventually, improve the stability of the plane in this range. 

3. STABILITY ANALYSIS. 

Consider linear time -invariant system 


x=Ax ( 4 ) 

where xeR n is a state vector. Then, assuming that the system is 
asymptotically stable one can define the following notions, [2]. 

Definition 1 . 

A connected set in the system parameters -space (parameters of 

matrix A) is a robust time invariant stable (RTIS) set for system 
(4) iff AeQ-j. and every time- invariant system 


x=<4x 


(5). 


is asymptotically stable for 

a 
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Definition 2. 


A connected set £2^. in the system parameters -space is a robust t: 
varying stable (RTVS) set for system (4) iff Ae£2y and every tii 
varying system (5) is asymptotically stable for 4Ze£2y. 


Then, consider four linear models instead of (2), respectively: 


-0.072815870 

for 

0° 

< 

a 

< 

14.74 

-0.048161261 

for 

14.74° 

< 

a 

< 

17.40 

-0.041453149 

for 

17.40° 

< 

a 

< 

18.87 

-0.016633734 

for 

18.87° 

< 

a 

< 

28 . 00 


It is easy to find that all models are asymptotically stable, 
are, however, interested in the set of (k^k^) such that all 
linear closed loop systems will be stable with the follow 
feedback 

u=Kx, K=[k 1 k£] 

An appropriate region Oj. can be easily calculated based 
algorithm 2 proposed in [2] . This is, however, only the sec 
order system and one can simply obtain analytical formulas for 
RTIS region in this case. The characteristic polynomial for 
2nd order system has the following form 

s^+as+b=0 

It is known that all roots of this polynomial are in the left b 
plane, i.e. a system is asymptotically stable stable, iff a>0 
b>0 . Based on this, the RTIS region for 'stable' feedback ga 
was calculated. The region, is presented on fig. 1, a dashed 1 
represents RTIS region for model PI, 0°sa<14 . 74° , a dotted moc 
P2 for 14 .74°<ail7 .4°, a dash-dotted model P3 for 17 . 40°<asl8 . 8 
and a continuous line model P4 for 18 . 87°<a^28° . It is easy to : 
that the system without feedback, i.e. k^=k 2 = 0 , sign + on 
plane (^^ 2 ), is very close to the stability region boundary, 
can improve stability assuming appropriate K from n I (k^,k 2 ). 


Next, RTVS sets f2y were calculated for these models, according to 
the algorithm given in [2], for uncertain parameters a^ and &21 
in A. They are presented on fig. 2. All four models are inside the 
RTVS region calculated for the model P4 . Moreover, since all time 
varying (nonlinear) a^ = 9.168c (a) is smaller than nominal values 
used in linear models it means that the whole nonlinear system (1) 
is asymptotically stable for 0°ia^28°. However, there is a very 
small upper bound for a^ in this model, namely 

+Aa^ < 0.0227 

This can cause that with small system uncertainty the system can 
be unstable. The vertexes of R TVS quadrilateral on the plane 
(Aa^,Aa 2 ^) are as follows: 

V VQ = { (-176.5,0), (0.0226,0), (0,-0.2274) (0,0.2944) } 

In order to improve system stability feedback gain matrix K was 
chosen from Q . Intuitively, it seems that a good gain is 'a small 
one - a high gain can result in a lack of system controllability 
because of saturation of the control input, and such that the 
stability margin with respect to K will be rather large. 

Thus, the good choice seems to be K^=[0 0.2]. For this gain one 
obtains significant improvement of RTVS set. This set is shown on 
fig. 3. In this case also all models are inside the calculated 
for the model P4 . However, an upper bound for a^ is more than 8 
times greater: 

+Aa^ < 0.1835 

Also range for uncertain parameter a 2 ^ is almost 6 times larger. 
Indeed, the vertexes of RTVS quadrilateral in this case on the 
plane ( Aa^ , Aa 2 -j_ ) are as follows 

V = { (-25049,0), (0.1835,0), (0,-3.547), (0,3.211) } 
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Then, it was considered feedback gain K2=[0.2 0]. This ga 

however, seems to be worse situated in the RTIS set than 

considering the stability region with respect to K. Neverthele 
also in this case one obtains improvement of robust stability 
the closed loop system. The appropriate RTVS set is shown 
fig. 4. In this case an upper bound for a^ is as follows 

+A a n < 0.0613 

Similarly, range for perturbation in a 21 is larger than for K 
The vertexes of RTVS set are as follows 

V v2 = { (-65.15,0), (0.0613,0), (0,-0.5770), (0,1.3641) } 

It should be noted that all RTVS sets were calculated uni 
assumption Q= I in algorithm 3 [2], 

From the above analysis follows that relatively small sta 
linear feedback gain K=[0 0.2] significantly improves stabil 

of the system. It should be emphsized that every nonlinear/ti: 
varying system (1) with a 11 and a 21 from the obtained RTVS set 
will be asymptotically stable. This way we have designed a robu 
stable nonlinear closed-loop system. 

4. CONCLUDING REMARKS. 

A robust-stable nonlinear controlsystem has been designed. It 
shown that small linear static feedback gain can significar 
improve stability of the airplane. The feedback gain seems to 
so small that it should not constrain control signal during pi 
maneuvering. This should also results in better a controllabil 
of the plane. 

/ 

The stability analysis and feedback gain synthesis were done us 
methods designed for linear systems [ 2 ] . This, approach can 
also used for more complicated nonlinear systems. For instar 
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assuming as a base model for the airplane, the linear 9th order 
model given in [1,9]. This model is unstable, but, as it was shown 
in [2], one can deal also with unstable models using the same 
approach. 

Presented results also show the power of the approach proposed in 

[ 2 ] . 
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